!     Copyright (C) 2010 Benjamin Piaud
!
!     LIMBES is free software; you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation; either version 3, or (at your option)
!     any later version.

!     LIMBES is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.

!     You should have received a copy of the GNU General Public License
!     along with LIMBES; if not, see <http://www.gnu.org/licenses/>

!     LIMBES is a recursive acronyme for 
!     LIMBES Is "Maybe" a Boltzmann Equation Solver.
!     LIMBES is based on algorithm between discrete velocities model family
!     and Lattice Boltzmann Methods

program LIMBES_main
  use LIMBES_mod_var
  use LIMBES_mod_fluid_config
  implicit none
  integer :: nb_periodes_initial,nb_periodes_final,nb_periodes_sec,nb_periodes_max,nb_periodes      
  real :: temps_elapsed
  integer*8 :: ni,ntot,time_index,itime,t_int,vi,nv
  real*8 :: dt

  !  time control init
  CALL SYSTEM_CLOCK(COUNT_RATE=nb_periodes_sec, COUNT_MAX=nb_periodes_max)
  CALL SYSTEM_CLOCK(COUNT=nb_periodes_initial)

  !     fetching simulation parameters and fluid properties
  call LIMBES_init()
  dt=SIMUL_dt       
  ntot=SIMUL_cell_number
  nv=LIMBES_nv
  itime=0
  time_index=0
  t_int=0

  !$omp parallel do
  do ni=1,ntot
     call LIMBES_macro(ni)
     call LIMBES_calc_feq(LIMBES_rho(ni),LIMBES_ux(ni),LIMBES_uy(ni),LIMBES_RT(ni),LIMBES_feq(ni,:))
  enddo

  call user_action(itime)

  !     boundary condition
  !$omp parallel do
  do ni=1,SIMUL_1wall_cell_number
     call LIMBES_boundary_1wall(ni)
  enddo
  !!!!!!!!$omp parallel do
  do ni=1,SIMUL_2wall_cell_number
     call LIMBES_boundary_2wall(ni)
  enddo
  !     saving initial macroscopic fields
  call LIMBES_write_macro(itime,time_index)
  call LIMBES_check_conserving(itime,time_index)
  call post_process(itime,time_index)

  !------------------------------------------------------------------
  !     main loop on time
  do itime=1,SIMUL_step_number
     !     loop on space to compute v.grad(f) and collision term
     !$omp parallel do
     do ni=1,ntot
        call LIMBES_transport(ni) ! stream step 
        call LIMBES_collision(ni) !collision step
     enddo

     !     boundary condition
     !$omp parallel do
     do ni=1,SIMUL_1wall_cell_number
        call LIMBES_boundary_1wall(ni)
     enddo
     !$omp parallel do
     do ni=1,SIMUL_2wall_cell_number
        call LIMBES_boundary_2wall(ni)
     enddo

     ! evolution equation: df/dt = - v.grad f + Coll
     LIMBES_f(:,:)=LIMBES_f(:,:) - LIMBES_vgradf(:,:)*dt + LIMBES_coll(:,:)*dt  

     !$omp parallel do
     do ni=1,ntot
        call LIMBES_macro(ni) ! update macroscopic fields 
        call LIMBES_calc_feq(LIMBES_rho(ni),LIMBES_ux(ni),LIMBES_uy(ni),LIMBES_RT(ni),LIMBES_feq(ni,:))
     enddo

     call user_action(itime)

     !     saving of results sometimes
     if (t_int.GE.SIMUL_step_number/SIMUL_save_number) then
        t_int=0
        time_index=time_index+1
        write(*,*) int(100*itime/(SIMUL_step_number))+1,'%'
        call LIMBES_write_macro(itime,time_index)
        call LIMBES_save_f()
        call LIMBES_check_conserving(itime,time_index)
        call post_process(itime,time_index)
     end if
     t_int=t_int+1

  enddo
  !     end of main loop on time
  !------------------------------------------------------------------

  !     last macroscopic fields saving
  call LIMBES_write_macro(SIMUL_step_number,time_index+1)
  call LIMBES_save_f()
  call LIMBES_check_conserving(SIMUL_step_number,time_index+1)  
  call post_process(SIMUL_step_number,time_index+1)  

  ! elapsed time
  CALL SYSTEM_CLOCK(COUNT=nb_periodes_final)
  nb_periodes = nb_periodes_final - nb_periodes_initial
  IF (nb_periodes_final < nb_periodes_initial) &
       nb_periodes = nb_periodes + nb_periodes_max
  temps_elapsed   = FLOAT (nb_periodes) / nb_periodes_sec
  print*,'elapsed time:', temps_elapsed

  ! free memory
  call free_LIMBES_var()
  call free_LIMBES_fluid_config()


end program LIMBES_main

!----------------------------------------------------------------
subroutine post_process(itime,time_index)
  use LIMBES_mod_var
  use LIMBES_mod_fluid_config
  implicit none
  integer*8 :: itime,time_index
!!!!! variables defined by the user !!!!
  integer*8 :: xi,yi,ni,vi,nwi,nwall
  real*8 :: w_f,Jtot,vx,vy
  
  Jtot=0.0D+0
  nwall=0  
  do nwi=1,SIMUL_1wall_cell_number         
     ni=SIMUL_1wall_cell(nwi)
     xi=LIMBES_cell_xi(ni)
     yi=LIMBES_cell_yi(ni)
     if (xi.EQ.496.AND.yi.LE.201) then
        nwall=nwall+1
        do vi=1,LIMBES_nv
           vx=LIMBES_vx(vi)
           vy=LIMBES_vy(vi)
           w_f=LIMBES_w(vi)*LIMBES_f_1wall(nwi,vi)
           Jtot=Jtot + 0.5D+0*vx*(vy*vy + vx*vx)*w_f         
        enddo
     endif
  enddo

  open(200,file='./data/J.dat')
  do ni=0,time_index-1
     read(200,*)
  enddo
  write(200,*) dble(itime)*SIMUL_dt,Jtot,Jtot/dble(nwall)
  close(200)
  
  
  open(200,file='./data/profil.dat')
  yi=250
  do xi=496,505
     ni=LIMBES_cell_ni(xi,yi)
     write(200,*) xi,LIMBES_ux(ni),LIMBES_uy(ni),LIMBES_RT(ni)/FLUID_R
  enddo
  close(200)



end subroutine post_process

!----------------------------------------------------------------

subroutine user_action(itime)
  use LIMBES_mod_var
  use LIMBES_mod_fluid_config
  implicit none
  integer*8 :: itime
!!!!! variables defined by the user !!!!
 

end subroutine user_action

